Magnet hypothesis: plant-pollinator interactions

Purpose: A test of the magnet hypothesis was examined in Mojave National Preserve by Ally Ruttan.

Hypothesis: Floral resource island created by shrubs and the associated beneficiary annual plants will positively and non-additively influence pollinator visitation rates.

Predictions:
(1) The frequency of pollinator visitations to annuals is greater under shrubs than in the paired-open microsites.
(2) Annual plants under flowering entomophilous shrubs (Larrea tridentata) will have a higher frequency of pollinator visitations than annual plants under anemophilous shrubs (Ambrosia dumosa) because of higher concentrations of floral resources for pollinators.
(3) Shrubs with annuals in their understory will have a higher frequency of pollinator visitations than shrubs without annuals due to increased concentrations of floral resources for pollinators.
(4) Sites with both shrubs and annuals will have the highest frequency of pollinator visitations to both the shrubs and the annuals.

An interesting corollary is that there are appropriate floral resources for desert pollinators, that they discriminate, and that entomophilous and anemophilous shrubs facilitate flowering similarly.

Data wrangling

#libraries
library(tidyverse)
library(DT)
library(lubridate)

#meta-data
meta <- read_csv("data/meta-data.csv")
datatable(meta)
#data
data.2015 <- read_csv("data/MNP.2015.csv")
data.2016 <- read_csv("data/MNP.2016.csv")

#merge
data <- rbind.data.frame(data.2015, data.2016)
data <- data %>% rename(net.treatment = treatment) %>% na.omit(data) 
data$year <- as.character(data$year)
data$rep <- as.character(data$rep)

#tidy data to expand treatment column (current structure is a mix of three factors)
#microsite
data <- data %>% mutate(microsite = ifelse(net.treatment %in% c("SA", "SAA", "SX"), "Larrea", ifelse(net.treatment %in% c("OA"), "open", ifelse(net.treatment %in% c("AMB"), "Ambrosia", NA))))
length(unique(data$microsite))
## [1] 3
#annuals
data <- data %>% mutate(annuals = ifelse(net.treatment %in% c("SA", "SAA"), "annuals", ifelse(net.treatment %in% c("OA"), "annuals", ifelse(net.treatment %in% c("AMB"), "annuals", "none"))))
length(unique(data$annuals))
## [1] 2
#flowering shrub
data <- data %>% mutate(flowering.shrub = ifelse(net.treatment %in% c("SAA", "SX"), "flowering shrub", ifelse(net.treatment %in% c("AMB", "SA", "OA"), "no shrub flowers", "NA")))
length(unique(data$annuals))
## [1] 2
#visitation duration needs to be weighted by duration (i.e. total recording time to address sampling effort)
data <- data %>% mutate(weighted.visitation.duration = as.duration(visitation.duration)/as.duration(total.duration))

#frequency counts in a separate dataframe (weighted by duration of recording)
frequency <- data %>% group_by(year, day, net.treatment, rep, microsite, annuals, flowering.shrub) %>% summarise(net.time = sum(total.duration), mean.visitation.duration = mean(weighted.visitation.duration), mean.floral.density = mean(floral.density), count = n()) 

frequency$net.time <- as.numeric(frequency$net.time) #converts to total seconds
frequency$rate <- as.numeric(frequency$count)/frequency$net.time #weight frequency by net time
datatable(frequency)

Data visualization

#higher-order treatment patterns in frequency####
ggplot(frequency, aes(net.treatment, rate)) + geom_boxplot() + facet_wrap(~year)

ggplot(frequency, aes(microsite, rate)) + geom_boxplot() + facet_wrap(~year*annuals)

ggplot(frequency, aes(microsite, rate)) + geom_boxplot() + facet_wrap(~year*annuals*flowering.shrub)

#relationships with sampling effort
ggplot(frequency, aes(net.time, count, color = year)) + geom_point()

ggplot(frequency, aes(net.time, count, color = year)) + geom_point() + facet_wrap(~microsite)

#floral density
ggplot(frequency, aes(mean.floral.density, rate, color = year)) + geom_point()

ggplot(frequency, aes(mean.floral.density, rate, color = microsite)) + geom_point() + facet_wrap(~year)

ggplot(frequency, aes(mean.floral.density, rate, color = microsite)) + geom_point() + facet_wrap(~year)

#mean visitation duration is also really important because you have more visits or longer visits####
ggplot(frequency, aes(net.treatment, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year)

ggplot(frequency, aes(microsite, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year*annuals)

ggplot(frequency, aes(microsite, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year*annuals*flowering.shrub)

#relationships with sampling effort
ggplot(frequency, aes(net.time, mean.visitation.duration, color = year)) + geom_point()

ggplot(frequency, aes(net.time, mean.visitation.duration, color = year)) + geom_point() + facet_wrap(~microsite)

#floral density
ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = year)) + geom_point()

ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = microsite)) + geom_point() + facet_wrap(~year)

ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = microsite)) + geom_point() + facet_wrap(~year)

EDA

#test distributions and explore outliers
summary(frequency)
##      year               day            net.treatment     
##  Length:266         Length:266         Length:266        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##      rep             microsite           annuals         
##  Length:266         Length:266         Length:266        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##  flowering.shrub       net.time      mean.visitation.duration
##  Length:266         Min.   :   300   Min.   :0.0000000       
##  Class :character   1st Qu.:  4500   1st Qu.:0.0008142       
##  Mode  :character   Median : 23010   Median :0.0065703       
##                     Mean   : 64855   Mean   :0.0159156       
##                     3rd Qu.:100660   3rd Qu.:0.0124122       
##                     Max.   :542929   Max.   :0.4941065       
##  mean.floral.density     count             rate          
##  Min.   :  5.00      Min.   :  1.00   Min.   :0.0001774  
##  1st Qu.: 23.13      1st Qu.:  3.00   1st Qu.:0.0002130  
##  Median : 71.26      Median :  8.00   Median :0.0002780  
##  Mean   : 98.74      Mean   : 15.14   Mean   :0.0005540  
##  3rd Qu.:200.00      3rd Qu.: 22.00   3rd Qu.:0.0011111  
##  Max.   :200.00      Max.   :109.00   Max.   :0.0033333
#check orthogonality
freq.2015 <- frequency %>% filter(year == 2015)
freq.2016 <- frequency %>% filter(year == 2016)

#annual treatment
ggplot(freq.2015, aes(rate, fill = annuals)) + geom_histogram() + facet_wrap(~microsite)

#only larrea tested annuals and no-annuals in 2015

ggplot(freq.2016, aes(rate, fill = annuals)) + geom_histogram() + facet_wrap(~microsite)

#only larrea tested annuals and no-annuals in 2016

#flowering shrubs
ggplot(freq.2015, aes(rate, fill = flowering.shrub)) + geom_histogram() + facet_wrap(~microsite)

#ha, only larrea in 2015

ggplot(freq.2016, aes(rate, fill = flowering.shrub)) + geom_histogram() + facet_wrap(~microsite)

#ok, only larrea again but had ambrosia

#conclusion, separate years
require(fitdistrplus)
descdist(freq.2015$rate, boot = 1000)

## summary statistics
## ------
## min:  0.0001873882   max:  0.001337793 
## median:  0.0003267998 
## mean:  0.0004887686 
## estimated sd:  0.0003435598 
## estimated skewness:  1.110484 
## estimated kurtosis:  2.810992
descdist(freq.2016$rate, boot = 1000)

## summary statistics
## ------
## min:  0.0001773679   max:  0.003333333 
## median:  0.0002480947 
## mean:  0.0006144842 
## estimated sd:  0.000504536 
## estimated skewness:  1.331961 
## estimated kurtosis:  7.146892
detach("package:fitdistrplus", unload = TRUE)

Models

#GLM for count and weight by net.time (alt - use MASS and glm.nb)
#library(MASS) #need for glm.nb

#all codes
#2015 counts
m <- glm(count~net.treatment + mean.floral.density, family = "poisson", weight = net.time, data = freq.2015)
anova(m, test = "Chisq") 
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: count
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  127   47238655              
## net.treatment        3 13736996       124   33501659 < 2.2e-16 ***
## mean.floral.density  1 13630462       123   19871197 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(freq.2015, aes(count, fill = microsite, weight = net.time)) + geom_histogram(binwidth = 10) + scale_fill_brewer(palette = "Blues") #like it

#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")
## $lsmeans
##  net.treatment   lsmean           SE df asymp.LCL asymp.UCL
##  OA            2.045110 0.0004169528 NA  2.044293  2.045927
##  SA            2.614917 0.0005830291 NA  2.613774  2.616060
##  SAA           1.736211 0.0004255844 NA  1.735377  1.737045
##  SX            1.354487 0.0007356275 NA  1.353045  1.355929
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate           SE df  z.ratio p.value
##  OA - SA  -0.5698067 0.0007644926 NA -745.340  <.0001
##  OA - SAA  0.3088991 0.0002022540 NA 1527.283  <.0001
##  OA - SX   0.6906229 0.0007383659 NA  935.340  <.0001
##  SA - SAA  0.8787058 0.0007683106 NA 1143.686  <.0001
##  SA - SX   1.2604296 0.0009583930 NA 1315.149  <.0001
##  SAA - SX  0.3817238 0.0007455663 NA  511.992  <.0001
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
#2016 counts
m <- glm(count~net.treatment + mean.floral.density, family = "poisson", weight = net.time, data = freq.2016)
anova(m, test = "Chisq") 
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: count
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  137  174945050              
## net.treatment        4 15411680       133  159533369 < 2.2e-16 ***
## mean.floral.density  1  8872062       132  150661307 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(freq.2016, aes(count, fill = microsite, weight = net.time)) + geom_histogram(binwidth = 10) + scale_fill_brewer(palette = "Blues") #good plot

#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")
## $lsmeans
##  net.treatment     lsmean           SE df  asymp.LCL  asymp.UCL
##  AMB            5.2074833 0.0005064262 NA  5.2064907  5.2084759
##  OA             5.2481438 0.0004770340 NA  5.2472089  5.2490788
##  SA            -1.0445768 0.0023971369 NA -1.0492751 -1.0398785
##  SAA            5.3522565 0.0004641314 NA  5.3513468  5.3531661
##  SX            -0.7880837 0.0019263237 NA -0.7918592 -0.7843082
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate           SE df   z.ratio p.value
##  AMB - OA  -0.04066053 1.176853e-04 NA  -345.502  <.0001
##  AMB - SA   6.25206012 2.597272e-03 NA  2407.164  <.0001
##  AMB - SAA -0.14477316 1.130315e-04 NA -1280.821  <.0001
##  AMB - SX   5.99556701 2.170318e-03 NA  2762.529  <.0001
##  OA - SA    6.29272065 2.583621e-03 NA  2435.621  <.0001
##  OA - SAA  -0.10411263 9.853844e-05 NA -1056.569  <.0001
##  OA - SX    6.03622754 2.153962e-03 NA  2802.383  <.0001
##  SA - SAA  -6.39683328 2.578069e-03 NA -2481.250  <.0001
##  SA - SX   -0.25649310 2.889407e-03 NA   -88.770  <.0001
##  SAA - SX   6.14034017 2.147300e-03 NA  2859.563  <.0001
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 5 estimates
#repeat above for mean.visitation.duration
#2015
m <- glm(mean.visitation.duration~net.treatment + mean.floral.density, family = "gaussian", weight = net.time, data = freq.2015)
anova(m, test = "Chisq") 
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: mean.visitation.duration
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                  127     2609.2              
## net.treatment        3   414.97       124     2194.2 3.553e-05 ***
## mean.floral.density  1     0.45       123     2193.7    0.8738    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")
## $lsmeans
##  net.treatment      lsmean          SE df    asymp.LCL  asymp.UCL
##  OA            0.013242840 0.005495694 NA  0.002471478 0.02401420
##  SA            0.001437038 0.007523933 NA -0.013309600 0.01618368
##  SAA           0.027506302 0.005617119 NA  0.016496950 0.03851565
##  SX            0.002694322 0.007347442 NA -0.011706400 0.01709504
## 
## Results are given on the identity (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate          SE df z.ratio p.value
##  OA - SA   0.011805802 0.010630339 NA   1.111  0.6830
##  OA - SAA -0.014263461 0.004063898 NA  -3.510  0.0025
##  OA - SX   0.010548518 0.008933348 NA   1.181  0.6390
##  SA - SAA -0.026069263 0.010652352 NA  -2.447  0.0685
##  SA - SX  -0.001257284 0.010632340 NA  -0.118  0.9994
##  SAA - SX  0.024811979 0.009016733 NA   2.752  0.0302
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
#2016
m <- glm(mean.visitation.duration~net.treatment + mean.floral.density, family = "gaussian", weight = net.time, data = freq.2016)
anova(m, test = "Chisq") 
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: mean.visitation.duration
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                                  137      57288           
## net.treatment        4   4296.8       133      52991   0.0300 *
## mean.floral.density  1     42.5       132      52949   0.7448  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")
## $lsmeans
##  net.treatment      lsmean         SE df  asymp.LCL asymp.UCL
##  AMB           -0.01369199 0.07143076 NA -0.1536937 0.1263097
##  OA             0.02699146 0.06771043 NA -0.1057185 0.1597015
##  SA             0.05059323 0.13300537 NA -0.2100925 0.3112790
##  SAA           -0.01396078 0.06610784 NA -0.1435298 0.1156082
##  SX             0.06860790 0.12744942 NA -0.1811884 0.3184042
## 
## Results are given on the identity (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast       estimate         SE df z.ratio p.value
##  AMB - OA  -0.0406834474 0.01550018 NA  -2.625  0.0659
##  AMB - SA  -0.0642852197 0.19379195 NA  -0.332  0.9974
##  AMB - SAA  0.0002687908 0.01515308 NA   0.018  1.0000
##  AMB - SX  -0.0822998902 0.19002170 NA  -0.433  0.9927
##  OA - SA   -0.0236017723 0.19045530 NA  -0.124  0.9999
##  OA - SAA   0.0409522382 0.01380555 NA   2.966  0.0251
##  OA - SX   -0.0416164428 0.18661765 NA  -0.223  0.9995
##  SA - SAA   0.0645540105 0.18909457 NA   0.341  0.9971
##  SA - SX   -0.0180146705 0.10981315 NA  -0.164  0.9998
##  SAA - SX  -0.0825686810 0.18522873 NA  -0.446  0.9918
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
#treatments split out
#how to test?

Intrepretation

Net.treatment models
1. Using the simplest models with all treatment levels aggregated but split by years, larrea is the most important magnet for frequency of visits (weighting for net.time recorded and using floral density as a covariate). This is a reasonable test because it captures enough of the variation and address non-orthogonality levels.

  1. Mean visitation rate in 2015 significantly differed between larea with flowers and open and larea with flowers and annuals versus larrea without. This strongly suggests a double-magnet effect this season. In 2016, larrea flowering with annuals was also significantly greater than open microsites with annuals.

  2. Double-magnet supported for frequency too.

  3. Then check each prediction. Need to think over how to handle double-magnet H.